// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtIndexSelectFilter.C
// ************************************************************************* //

#include <avtIndexSelectFilter.h>

#include <avtExecutionManager.h>

#include <vtkCell.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCharArray.h>
#include <vtkDataSetReader.h>
#include <vtkDataSetRemoveGhostCells.h>
#include <vtkDataSetWriter.h>
#include <vtkMaskPoints.h>
#include <vtkIntArray.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnsignedIntArray.h>
#include <vtkUnstructuredGrid.h>
#include <vtkVisItExtractGrid.h>
#include <vtkVisItExtractRectilinearGrid.h>
#include <vtkVisItUtility.h>

#include <avtCallback.h>
#include <avtDatasetExaminer.h>
#include <avtLogicalSelection.h>
#include <avtSILNamespace.h>
#include <avtSILRestrictionTraverser.h>
#include <avtOriginatingSource.h>
#include <avtParallel.h>

#include <DebugStream.h>
#include <ImproperUseException.h>
#include <InvalidSetException.h>
#include <InvalidCategoryException.h>
#include <InvalidVariableException.h>

#ifdef PARALLEL
#include <mpi.h>
#endif

#include <string>
#include <vector>
using std::string;
using std::vector;

// ****************************************************************************
//  Method: avtIndexSelectFilter constructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Jun 5 09:09:11 PDT 2002
//
//  Modifications:
//
//    Mark C. Miller, Tue Sep 28 19:57:42 PDT 2004
//    Added data selection id
//
//    Kathleen Bonnell, Tue May 10 11:19:24 PDT 2005 
//    Use VisIt versions of vtkExtractGrid and vtkExtractRectilinearGrid, 
//    they have been modified to correctly handle cell data when VOI is
//    along max boundary. 
//
//    Kathleen Bonnell, Mon Jan 30 15:10:26 PST 2006 
//    Add vtkMaskPoints for a points filter. 
// 
//    Kathleen Bonnell, Thu Jun  7 14:37:49 PDT 2007 
//    Added groupCategory.
// 
//    Kathleen Bonnell, Thu Jun 21 16:31:59 PDT 2007 
//    Added amrLevel, amrMesh.
//
//    Kathleen Biagas, Wed Jan 11 13:58:31 PST 2012
//    Turn on SingleCellPerVertex for vtkMaskPoints.
//
//    Kathleen Biagas, Tue Jun 9 09:31:15 MST 2015
//    Added globalDims, lspace.
//
//    Alister  Maguire, Wed Oct 26 11:34:15 PDT 2016
//    Removed curvilinearFilter, rectilinearFilter, 
//    and pointsFilter for thread safety. They now
//    exist as stack variables in ExecuteData. 
//
// ****************************************************************************

avtIndexSelectFilter::avtIndexSelectFilter() : lspace()
{
    haveIssuedWarning = false;
    selID             = -1;
    amrLevel          = -1;
    amrMesh           = false;
    groupCategory = false;
    for (int i = 0; i < 3; ++i)
    {
        globalDims[i] = INT_MAX;
        globalDims[i+3] = -1;
    }
}


// ****************************************************************************
//  Method: avtIndexSelectFilter destructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Jun 5 09:09:11 PDT 2002
//
//  Modifications:
//    Kathleen Bonnell, Mon Jan 30 15:10:26 PST 2006 
//    Delete vtkMaskPoints.
//
//    Alister  Maguire, Wed Oct 26 11:34:15 PDT 2016
//    Removed curvilinearFilter, rectilinearFilter, 
//    and pointsFilter for thread safety. They now
//    exist as stack variables in ExecuteData. 
//
// ****************************************************************************

avtIndexSelectFilter::~avtIndexSelectFilter()
{

}


// ****************************************************************************
//  Method:  avtIndexSelectFilter::Create
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Jun 5 09:09:11 PDT 2002
//
// ****************************************************************************

avtFilter *
avtIndexSelectFilter::Create()
{
    return new avtIndexSelectFilter();
}


// ****************************************************************************
//  Method:      avtIndexSelectFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Jun 5 09:09:11 PDT 2002
//
// ****************************************************************************

void
avtIndexSelectFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const IndexSelectAttributes*)a;
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::PrepareFilters
//
//  Purpose:
//      Prepares the curvilinear and rectilinear filters for the current
//      dataset.  This can take into account block indices.
//
//  Programmer: Hank Childs
//  Creation:   June 25, 2002
//
//  Modifications:
//    Kathleen Bonnell, Wed Sep  8 09:36:30 PDT 2004
//    Always Set IncludeBoundary to true for filters, so they can handle
//    modulo prolbems (eg sample rate of 3, but dimension is 10).
//
//    Kathleen Bonnell, Wed Jul 20 11:39:34 PDT 2005 
//    Don't subtract 1 from the groupIndices. 
//
//    Kathleen Bonnell, Mon Jan 30 15:10:26 PST 2006 
//    Setup vtkMaskPoints.
//
//    Kathleen Bonnell, Thu Jun  7 14:37:49 PDT 2007 
//    Use groupCategory insead of WhichData == OneGroup.
// 
//    Kathleen Bonnell, Thu Jun 21 16:31:59 PDT 2007 
//    Added int* amri arg.  Read from this arg if AMR mesh in order to
//    determine how to correctly compute min/max.
//
//    Alister Maguire, Mon Oct 24 12:25:39 PDT 2016
//    Added vtkVisItExtgracGrid, vtkVisItExtgractRectilinearGrid, 
//    and vtkMaskPoints arguments for thread safety.
//
// ****************************************************************************

void
avtIndexSelectFilter::PrepareFilters(int groupIndices[3], int *amri, 
                                     vtkVisItExtractGrid *curvilinearFilter,
                                     vtkVisItExtractRectilinearGrid *rectilinearFilter, 
                                     vtkMaskPoints *pointsFilter)
{
    //
    // Only adjust the base index if we are index selecting by group.
    //
    int bi[3];
    if (groupCategory)
    {
        bi[0] = groupIndices[0];
        bi[1] = groupIndices[1];
        bi[2] = groupIndices[2];
    }
    else
    {
        bi[0] = 0;
        bi[1] = 0;
        bi[2] = 0;
    }

    int voi[6];

    int minmax[6] = {atts.GetXMin(), atts.GetXMax(), atts.GetYMin(),
                     atts.GetYMax(), atts.GetZMin(), atts.GetZMax()};

    if (amrMesh && groupCategory)
    {
        int OP = amri[3];
        for (int i = 0; i < 3; i++)
        {
            if (amri[i] > 1)
            {
                if (OP == 0)
                {
                    // min
                    minmax[i*2] = minmax[i*2]*amri[i]; 
                    // max
                    minmax[i*2+1] = minmax[i*2+1]*amri[i]; 
                }
                else 
                {
                    // min
                    minmax[i*2] = minmax[i*2]/amri[i]; 
                    // max
                    if (minmax[i*2+1] != -1)
                    minmax[i*2+1] = minmax[i*2+1]/amri[i] + 1; 
                }
            }
        }
    }
  
    voi[0] = (minmax[0] - bi[0] < 0 ? 0 : minmax[0] - bi[0]);
    voi[1] = (minmax[1] < 0 ? 1000000   : minmax[1] - bi[0]);
    if (atts.GetDim() == IndexSelectAttributes::TwoD ||
        atts.GetDim() == IndexSelectAttributes::ThreeD)
    {
        voi[2] = (minmax[2] - bi[1] < 0 ? 0 : minmax[2] - bi[1]);
        voi[3] = (minmax[3] < 0 ? 1000000   : minmax[3] - bi[1]);
    }
    else
    {
        voi[2] = 0;
        voi[3] = 1000000;
    }
    if (atts.GetDim() == IndexSelectAttributes::ThreeD)
    {
        voi[4] = (minmax[4] - bi[2] < 0 ? 0 : minmax[4] - bi[2]);
        voi[5] = (minmax[5] < 0 ? 1000000   : minmax[5] - bi[2]);
    }
    else
    {
        voi[4] = 0;
        voi[5] = 1000000;
    }

    curvilinearFilter->SetVOI(voi);
    rectilinearFilter->SetVOI(voi);
    int sampleRate[3];
    sampleRate[0] = atts.GetXIncr();
    sampleRate[1] = atts.GetYIncr();
    sampleRate[2] = atts.GetZIncr();
    curvilinearFilter->SetSampleRate(sampleRate);
    rectilinearFilter->SetSampleRate(sampleRate);

    pointsFilter->SetOnRatio(sampleRate[0]);
    pointsFilter->SetOffset(voi[0]);
    if (voi[1] != 1000000)
    {
        int maxpts = (voi[1] - voi[0]) / sampleRate[0]; 
        pointsFilter->SetMaximumNumberOfPoints(maxpts);
    }
    else
    {
        pointsFilter->SetMaximumNumberOfPoints(VTK_INT_MAX);
    }
    curvilinearFilter->SetIncludeBoundary(1);
    rectilinearFilter->SetIncludeBoundary(1);
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtIndexSelectFilter with the given
//      parameters would result in an equivalent avtIndexSelectFilter.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Jun 5 09:09:11 PDT 2002
//
// ****************************************************************************

bool
avtIndexSelectFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(IndexSelectAttributes*)a);
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the IndexSelect filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Wed Jun 5 09:09:11 PDT 2002
//
//  Modifications:
//
//    Hank Childs, Wed Jun 19 09:47:01 PDT 2002
//    Fix stupid crash when applying to unstructured grids.
//
//    Hank Childs, Tue Jun 25 18:22:35 PDT 2002
//    Use group indices when calculating indices.
//
//    Hank Childs, Wed Oct  2 10:13:26 PDT 2002
//    Added support for meshes that have already been broken up.
//
//    Hank Childs, Mon Dec  9 08:53:15 PST 2002
//    Account for potentially changing format of avtOriginalCellNumbers.
//
//    Kathleen Bonnell, Fri Dec 13 16:41:12 PST 2002   
//    Use NewInstance instead of MakeObject in order to match vtk's new api. 
//    
//    Jeremy Meredith, Fri Jan 30 17:45:47 PST 2004
//    Added code to preserve dataset type if the input was polydata.
//
//    Hank Childs, Fri Aug 27 15:25:22 PDT 2004
//    Rename ghost data array.
//
//    Mark C. Miller, Tue Sep 28 19:57:42 PDT 2004
//    Added code to bypass the operator if the selection is applied by a
//    plugin
//
//    Kathleen Bonnell, Fri Feb 18 09:41:16 PST 2005 
//    Account for the fact the vtkExtractGrid and vtkRectilinearExtractGrid
//    may return 'empty' datasets, so we want to return a NULL dataset. 
//
//    Kathleen Bonnell, Fri Aug 19 15:45:46 PDT 2005
//    Remove ghost nodes array, ensure ghost zones are skipped in
//    'zones not preserved' case.   Removed retrieval and use of
//    'avtRealDims'.  Its retrieval was incorrect so it was never ever used.
//
//    Kathleen Bonnell, Mon Jan 30 15:10:26 PST 2006 
//    Use vtkMaskPoints for point meshes.
//
//    Kathleen Bonnell, Wed May 17 10:46:58 PDT 2006
//    Remove call to SetSource(NULL) as it now removes information necessary
//    for the dataset. 
//
//    Kathleen Bonnell, Thu Jun  7 14:37:49 PDT 2007 
//    Use groupCategory insead of WhichData == OneGroup.
// 
//    Kathleen Bonnell, Thu Jun 21 16:31:59 PDT 2007 
//    If this is an AMR mesh, retrieve AMR indices and don't remove ghost zones.
//
//    Kathleen Bonnell, Tue Feb 22 11:44:48 PST 2011
//    Added code to convert vtkPolyVertex cells returned from vtkMaskPoints
//    to vtkVertex cells, as filters down the pipeline may not be able to
//    handle vtkPolyVertexCells.
//
//    Kathleen Biagas, Wed Jan 11 13:59:21 PST 2012
//    Remove hack that converts vtkPolyVertex cells to vtkVertex cells.
//
//    Eric Brugger, Mon Jul 28 15:33:52 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Eric Brugger, Fri Sep 26 08:48:14 PDT 2014
//    I modified the routine to return a NULL in the case where it previously
//    returned an avtDataRepresentation with a NULL vtkDataSet.
//
//    Kathleen Biagas, Tue Jun 9 09:34:17 MST 2015
//    Move Replicate (wrap) to PostExecute, collect globalDims here instead.
//
//    Alister Maguire, Tue Nov  8 12:44:55 PST 2016
//    Changed curvilinearFilter, rectilinearFilter, and pointsFilter to stack
//    variables for thread safety; added amrMesh for thread safety; added
//    mutex locks where necessary.  
//
// ****************************************************************************

avtDataRepresentation *
avtIndexSelectFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    
    vtkVisItExtractGrid            *curvilinearFilter = 
                                     vtkVisItExtractGrid::New();
    vtkVisItExtractRectilinearGrid *rectilinearFilter = 
                                     vtkVisItExtractRectilinearGrid::New();
    vtkMaskPoints                  *pointsFilter      = vtkMaskPoints::New();
    pointsFilter->GenerateVerticesOn();
    pointsFilter->RandomModeOff();
    pointsFilter->SingleVertexPerCellOn();

    //
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    vtkDataSet *out_ds = NULL;

    int topoDim = GetInput()->GetInfo().GetAttributes().
                  GetTopologicalDimension();

    bool amrMesh = (GetInput()->GetInfo().GetAttributes().GetMeshType() == 
                     AVT_AMR_MESH);   

    vtkIntArray *bi_arr = vtkIntArray::SafeDownCast(
        in_ds->GetFieldData()->GetArray("base_index"));
    //
    // If the selection this filter exists to create has already been handled,
    // then we can skip execution
    //
    if (GetInput()->GetInfo().GetAttributes().GetSelectionApplied(selID)) 
    {
        debug1 << "Bypassing IndexSelect operator because database plugin "
                  "claims to have applied the selection already" << endl;

        out_ds = in_ds;
    }
    else if (GetInput()->GetInfo().GetValidity().GetZonesPreserved())
    {
        //
        // We have the normal case -- a structured mesh that we are going to
        // index select.
        //
        vtkDataSet *ds = NULL;

        //
        // All of our indices are incorrect if we leave the ghost zones in -- 
        // we can also have some weird phenomenon where boundaries are missing
        // between blocks.
        //
        vtkDataSetRemoveGhostCells *removeGhostCells = NULL;
        if (!amrMesh &&  
            in_ds->GetCellData()->GetArray("avtGhostZones"))
        {
            removeGhostCells  = vtkDataSetRemoveGhostCells::New();
            removeGhostCells->SetInputData(in_ds);

            //
            // There is something buggy about the extents when this filter is 
            // used for repeated executions.  Just force the execution now.
            //
            removeGhostCells->Update();
            ds = removeGhostCells->GetOutput()->NewInstance();
            ds->ShallowCopy(removeGhostCells->GetOutput());
        }
        else
        {
            ds = in_ds;
        }
        in_ds->GetPointData()->RemoveArray("avtGhostNodes");
    
        //
        // The indices should reflect the "base_index"'s, so dummy one up if
        // we don't have one.
        //
        int ind[3] = { 0, 0, 0 };
        if (bi_arr != NULL)
        {
            ind[0] = bi_arr->GetValue(0);
            ind[1] = bi_arr->GetValue(1);
            ind[2] = bi_arr->GetValue(2);
        }

        int *amri = NULL;
        if (amrMesh && groupCategory)
        {
            vtkIntArray *amrdims = (vtkIntArray *)
                in_ds->GetFieldData()->GetArray("avtAMRDimensions");
            if (amrdims == NULL)
            {
                if (!haveIssuedWarning)
                {
                    avtCallback::IssueWarning("An internal error occurred and "
                            "the index select operator was not applied.");
                    VisitMutexLock("avtIndexSelectFilter::ThreadSafeExecuteData");
                    haveIssuedWarning = true;
                    VisitMutexUnlock("avtIndexSelectFilter::ThreadSafeExecuteData");
                }
                return in_dr;
            }
            else
            {
                amri = (int*)amrdims->GetVoidPointer(0);
            }
        }
        PrepareFilters(ind, amri, curvilinearFilter, rectilinearFilter, 
                       pointsFilter);
    
        vtkDataSet *rv = NULL;
        int dstype = ds->GetDataObjectType();
        if (dstype == VTK_STRUCTURED_GRID)
        {
            curvilinearFilter->SetInputData(ds);
            curvilinearFilter->Update();
            rv = curvilinearFilter->GetOutput();
        }
        else if (dstype == VTK_RECTILINEAR_GRID)
        {
            rectilinearFilter->SetInputData(ds);
            rectilinearFilter->Update();
            rv = rectilinearFilter->GetOutput();
        }
        else if (topoDim == 0 &&
                 (dstype == VTK_POLY_DATA || dstype == VTK_UNSTRUCTURED_GRID))
        {
            pointsFilter->SetInputData(ds);
            pointsFilter->Update();
            rv = pointsFilter->GetOutput();
        }
        else
        {
            if (!haveIssuedWarning)
            {
                avtCallback::IssueWarning("The index select operator was "
                            "applied to a non-structured mesh.  It is not "
                            "being applied.");
                VisitMutexLock("avtIndexSelectFilter::ThreadSafeExecuteData");
                haveIssuedWarning = true; 
                VisitMutexUnlock("avtIndexSelectFilter::ThreadSafeExecuteData");
            }
            return in_dr;
        }

        if (removeGhostCells != NULL)
        {
            removeGhostCells->Delete();
        }

        if (rv->GetNumberOfPoints() <= 0 || rv->GetNumberOfCells() <= 0)
        {
            return NULL;
        }

        out_ds = (vtkDataSet *) rv->NewInstance();
        out_ds->ShallowCopy(rv);
    }
    else
    {
        //
        // The dataset has been changed before it got here -- most likely it
        // was material selected.  We should have passed enough clues
        // downstream to figure out what happened.
        //
        vtkUnsignedIntArray *origZones = (vtkUnsignedIntArray *)
                                         in_ds->GetCellData()->
                                         GetArray("avtOriginalCellNumbers"); 
        if (origZones == NULL)
        {
            if (!haveIssuedWarning)
            {
                avtCallback::IssueWarning("An internal error occurred and the "
                                          "index select operator was not "
                                          "applied.");
                VisitMutexLock("avtIndexSelectFilter::ThreadSafeExecuteData");
                haveIssuedWarning = true;
                VisitMutexUnlock("avtIndexSelectFilter::ThreadSafeExecuteData");
            }
            return in_dr;
        }

        vtkUnsignedIntArray *dims = (vtkUnsignedIntArray *)
                                            in_ds->GetFieldData()->GetArray(
                                            "avtOriginalStructuredDimensions"); 
        if (dims == NULL)
        {
            if (!haveIssuedWarning)
            {
                avtCallback::IssueWarning("An internal error occurred and the "
                                          "index select operator was not "
                                          "applied.");
                VisitMutexLock("avtIndexSelectFilter::ThreadSafeExecuteData");
                haveIssuedWarning = true; 
                VisitMutexUnlock("avtIndexSelectFilter::ThreadSafeExecuteData");
            }
            return in_dr;
        }
        int d[3];
        d[0] = dims->GetValue(0);
        d[1] = dims->GetValue(1);
        d[2] = dims->GetValue(2);

        int base[3] = { 0, 0, 0 };
        if (groupCategory)
        {
            if (bi_arr != NULL)
            {
                base[0] = bi_arr->GetValue(0);
                base[1] = bi_arr->GetValue(1);
                base[2] = bi_arr->GetValue(2);
            }
        }
        
        vtkDataArray *ghosts = in_ds->GetCellData()->GetArray("avtGhostZones");
        unsigned char *gz = NULL;
        if (ghosts)
            gz = (unsigned char *)ghosts->GetVoidPointer(0);

        //
        // We should have everything lined up now -- we know what the original
        // indexing of the structured mesh was and what the base index is.
        //
        vtkUnstructuredGrid *out_ug = vtkUnstructuredGrid::New();
        vtkPoints *p1 = vtkVisItUtility::GetPoints(in_ds);
        out_ug->SetPoints(p1);
        p1->Delete();
        out_ug->GetPointData()->PassData(in_ds->GetPointData());
        vtkCellData *out_cd = out_ug->GetCellData();
        vtkCellData *in_cd  = in_ds->GetCellData();
        out_cd->CopyAllocate(in_cd);
        int ncells = in_ds->GetNumberOfCells();
        out_ug->Allocate(ncells);

        int out_cell = 0;
        int xmin = atts.GetXMin();
        int xmax = atts.GetXMax();
        xmax = (xmax < 0 ? 10000000 : xmax);
        int ymin = atts.GetYMin();
        int ymax = atts.GetYMax();
        ymax = (ymax < 0 ? 10000000 : ymax);
        int zmin = atts.GetZMin();
        int zmax = atts.GetZMax();
        zmax = (zmax < 0 ? 10000000 : zmax);
        int ncomps = origZones->GetNumberOfComponents();
        int comp = ncomps-1;
        for (int i = 0 ; i < ncells ; i++)
        {
            int cell_id = (int) origZones->GetComponent(i, comp);
            int x = cell_id % (d[0]-1);
            int y = (cell_id / (d[0]-1)) % (d[1]-1);
            int z = cell_id / ((d[0]-1)*(d[1]-1));
            x += base[0];
            y += base[1];
            z += base[2];
            if (x < xmin || x >= xmax)
            {
                continue;
            }
            if (atts.GetDim() == IndexSelectAttributes::TwoD ||
                atts.GetDim() == IndexSelectAttributes::ThreeD)
            {
                if (y < ymin || y >= ymax)
                    continue;
            }
            if (atts.GetDim() == IndexSelectAttributes::ThreeD)
            {
                if (z < zmin || z >= zmax)
                    continue;
            }
            // only allow AMR ghosts to be part of the output.
            if (gz && !(gz[i] == 0 || gz[i] == 8))
                continue;

            out_cd->CopyData(in_cd, i, out_cell++);
            vtkCell *cell = in_ds->GetCell(i);
            vtkIdList *list = cell->GetPointIds();
            out_ug->InsertNextCell(in_ds->GetCellType(i), list);
        }
        out_ds = out_ug;

        if (out_ds->GetNumberOfCells() <= 0)
        {
            out_ds->Delete();
            return NULL;
        }

        //
        // If we had poly data input, we want poly data output.  The VTK filter
        // only returns unstructured grids, so convert that now.
        //
        if (in_ds->GetDataObjectType() == VTK_POLY_DATA)
        {
            vtkUnstructuredGrid *ugrid = (vtkUnstructuredGrid *) out_ds;
            vtkPolyData *out_pd = vtkPolyData::New();
            out_pd->SetPoints(ugrid->GetPoints());
            out_pd->GetPointData()->ShallowCopy(ugrid->GetPointData());
            out_pd->GetCellData()->ShallowCopy(ugrid->GetCellData());
            int ncells = ugrid->GetNumberOfCells();
            out_pd->Allocate(ncells);
            for (int i = 0 ; i < ncells ; i++)
            {
                int celltype = ugrid->GetCellType(i);
                vtkIdType npts;
                const vtkIdType *pts;
                ugrid->GetCellPoints(i, npts, pts);
                out_pd->InsertNextCell(celltype, npts, pts);
            }
            out_ds->Delete();
            out_ds = out_pd;
        }
    }

    VisitMutexLock("avtIndexSelectFilter::ExecuteData");
    atLeastOneThreadSuccessfullyExecuted = true;
    VisitMutexUnlock("avtIndexSelectFilter::ExecuteData");

    if ((atts.GetXWrap() || atts.GetYWrap() || atts.GetZWrap()) &&
        topoDim > 0 &&
       (in_ds->GetDataObjectType() == VTK_STRUCTURED_GRID ||
        in_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID))
    {
        // In the multi-block case, this doesn't work if sil-selection
        // has been applied, because only the currently plotted domains
        // are used to calculate the maximums... is this okay?

        int lo[3] = {0, 0, 0};
        int hi[3] = {0, 0, 0};
        int dims[3] = {0, 0, 0};
        if (out_ds->GetDataObjectType() == VTK_STRUCTURED_GRID)
        {
            vtkStructuredGrid::SafeDownCast(out_ds)->GetDimensions(dims);
        }
        else if (out_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        {
            vtkRectilinearGrid::SafeDownCast(out_ds)->GetDimensions(dims);
        }
        if (bi_arr)
        {
            lo[0] = bi_arr->GetValue(0);
            lo[1] = bi_arr->GetValue(1);
            lo[2] = bi_arr->GetValue(2);
        }
        for (int i = 0; i < 3; ++i)
        {
            hi[i] = lo[i] + dims[i] -1;
            if (lo[i] < globalDims[i])
                globalDims[i] = lo[i];
            if (hi[i] > globalDims[i+3])
                globalDims[i+3] = hi[i];
        }
    }

    avtDataRepresentation *out_dr = new avtDataRepresentation(out_ds,
          in_dr->GetDomain(), in_dr->GetLabel());

    if (out_ds != in_ds)
        out_ds->Delete();

    if (curvilinearFilter != NULL)
    {
        curvilinearFilter->Delete();
        curvilinearFilter = NULL;
    }
    if (rectilinearFilter != NULL)
    {
        rectilinearFilter->Delete();
        rectilinearFilter = NULL;
    }
    if (pointsFilter != NULL)
    {
        pointsFilter->Delete();
        pointsFilter = NULL;
    }
    
    return out_dr;
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::ModifyContract
//
//  Purpose:
//      Restricts the SIL to the domains requested by the user.
//
//  Programmer: Hank Childs
//  Creation:   June 5, 2002
//
//  Modifications:
//
//    Hank Childs, Sun Jun 16 20:50:53 PDT 2002
//    Add support for non 0-origin blocks.
//
//    Hank Childs, Mon Sep 30 17:23:33 PDT 2002
//    Add support for index selecting after a destructive operation.
//
//    Hank Childs, Mon Dec  2 09:59:56 PST 2002
//    Account for changing interface for SIL restriction.
//
//    Hank Childs, Thu Aug 14 07:44:49 PDT 2003
//    Also request the structured indices if we are specifically told to do
//    interface reconstruction.
//
//    Mark C. Miller, Tue Sep 28 19:57:42 PDT 2004
//    Added code to build a data selection
//
//    Kathleen Bonnell, Tue Nov 16 16:13:08 PST 2004 
//    Gracefully handle domainIndex that is out-of-range, and issue warning. 
//    Also, use domainIndex when determining chunk for trav.GetMaterials, when
//    appropriate.
//
//    Kathleen Bonnell, Thu Jul 21 07:52:49 PDT 2005 
//    For OneGroup selection, determine domains belonging to that group,
//    and restrict by domains to ensure that sets that were previously
//    turned off remains so.  When determining if StructuredIndices are
//    required, ensure checking is done only on domains needed by this filter. 
//
//    Kathleen Bonnell, Thu Aug  4 15:47:59 PDT 2005 
//    Request original node numbers when required. 
//    
//    Jeremy Meredith, Wed Aug 24 13:37:07 PDT 2005
//    Added support for group origin.
//
//    Kathleen Bonnell, Thu Jun  7 14:37:49 PDT 2007 
//    Modified how SIL selection is done, based on new atts. 
// 
//    Kathleen Bonnell, Thu Jun 21 16:31:59 PDT 2007
//    If this is an AMR mesh and category is group, don't restrict the SIL
//    and also request AMR indices.
// 
//    Hank Childs, Fri Oct 26 16:40:57 PDT 2007
//    Correct some SIL handling.  Compact SIL attributes are only defined
//    over a subset of the nodes in the SIL graph; comparing them without
//    re-indexing ... which was being done previously ... led to incorrect
//    results.
//
//    Hank Childs, Wed Jan  9 16:10:33 PST 2008
//    Beef up logic to handle species selection.  We were turning those sets
//    off, which resulted in all zeroes.
//
//    Dave Bremer, Thu Jan 31 17:52:55 PST 2008
//    Small tweak to guard against a case in which the RealMapsOut are 
//    requested from an avtSILSet, but the set goes out of scope and its maps
//    out are deleted before this method is done using them.
//
//    Hank Childs, Mon Dec 14 16:55:10 PST 2009
//    Add support for new SIL interface.
//
// ****************************************************************************

avtContract_p
avtIndexSelectFilter::ModifyContract(avtContract_p spec)
{
    avtContract_p rv = new avtContract(spec);

    amrMesh = GetInput()->GetInfo().GetAttributes().GetMeshType() == AVT_AMR_MESH;
    bool skipSILRestriction = amrMesh && groupCategory;

    if (!atts.GetUseWholeCollection() && !skipSILRestriction) 
    {
        std::string category = atts.GetCategoryName();
        std::string subset = atts.GetSubsetName();
        avtSILRestriction_p silr = spec->GetDataRequest()->GetRestriction();
        avtSILRestriction_p old_values = new avtSILRestriction(silr);
        avtSILRestrictionTraverser trav(old_values);
        int collectionID = silr->GetCollectionIndex(category, silr->GetTopSet());
        int setID = silr->GetSetIndex(subset, collectionID);
        if (trav.UsesSetData(setID) == NoneUsed) 
        {
            EXCEPTION1(InvalidSetException, subset.c_str());
        }
        TRY
        {
            // If we've got species info, we need to maintain that.
            // So see which species are on.
            vector<int>  species;
            vector<bool> setState;
        
            int topset = silr->GetTopSet();
            avtSILSet_p pTopset = silr->GetSILSet(topset);

            const vector<int> &mapsOut = pTopset->GetRealMapsOut();
            avtSILCollection_p speciesColl = NULL;
            for (size_t i = 0 ; i < mapsOut.size() ; i++)
            {
                avtSILCollection_p coll = silr->GetSILCollection(mapsOut[i]);
                if (coll->GetRole() == SIL_SPECIES)
                {
                    speciesColl = coll;
                }
            }
            if (*speciesColl != NULL)
            {
                for (int i = 0 ; i < speciesColl->GetNumberOfSubsets() ; i++)
                    setState.push_back(trav.UsesData(speciesColl->GetSubset(i)));
            }
            // End logic for seeing which species is on.

            silr = rv->GetDataRequest()->GetRestriction();
            silr->TurnOffAll();
            silr->TurnOnSet(setID);
            // We've just turned on an entire set, but some parts
            // (materials) may have been turned off before, so ensure
            // that remains the case.
            int numSets = silr->GetNumSets();
            for (int i = 0; i < numSets ; i++)
            {
                if (setID == i)
                    continue;
                if (trav.UsesSetData(i) == NoneUsed)
                    silr->TurnOffSet(i);
            }

            // Turn sets back on if species are on.
            for (size_t i = 0 ; i < species.size() ; i++)
                if (setState[i])
                    silr->TurnOnSet(species[i]);
        }
        CATCH(InvalidVariableException)
        {
            // If for some reason the GetSetIndex fails.
            RETHROW;
        }
        ENDTRY
    }

    if (amrMesh && groupCategory && amrLevel != -1)
    {
        rv->GetDataRequest()->SetNeedAMRIndices(amrLevel);
    }

    if (!GetInput()->GetInfo().GetValidity().GetZonesPreserved())
    {
        rv->GetDataRequest()->SetNeedStructuredIndices(true);
    }
    else if (rv->GetDataRequest()->MustDoMaterialInterfaceReconstruction())
    {
        rv->GetDataRequest()->SetNeedStructuredIndices(true);
    }
    else 
    {
        bool needSI = false;
        avtSILRestriction_p silr =rv->GetDataRequest()->GetRestriction();
        avtSILRestrictionTraverser trav(silr);
        if (atts.GetUseWholeCollection())
        {
            intVector chunks;
            trav.GetDomainList(chunks);
            for (size_t i = 0; i < chunks.size(); i++)
            {
                bool hasMats = false;
                trav.GetMaterials(chunks[i], hasMats);
                needSI |= hasMats;
            }
        }
        else 
        {
            needSI = !trav.UsesAllMaterials();
        }
        if (needSI)
        {
            rv->GetDataRequest()->SetNeedStructuredIndices(true);
        }
    }

    //
    // Indicate this operator's data selection 
    //
    avtLogicalSelection *sel = new avtLogicalSelection;
    switch (atts.GetDim())
    {
        case IndexSelectAttributes::OneD:   sel->SetNDims(1); break;
        case IndexSelectAttributes::TwoD:   sel->SetNDims(2); break;
        case IndexSelectAttributes::ThreeD: sel->SetNDims(3); break;
    }
    int vec[3];
    vec[0] = atts.GetXMin();
    vec[1] = atts.GetYMin();
    vec[2] = atts.GetZMin();
    sel->SetStarts(vec);

    // IndexSelect is nodal based and is inclusive. As such, the user
    // can have the min == max. Which would give a slice from a volume.

    // avtLogicalSelection's stops are nodal and are inclusive
    // also, we need to deal with using '-1' to mean 'max'
    if (atts.GetXMax() == -1)
        vec[0] = -1;
    else
        vec[0] = atts.GetXMax();

    if (atts.GetYMax() == -1)
        vec[1] = -1;
    else
        vec[1] = atts.GetYMax();

    if (atts.GetZMax() == -1)
        vec[2] = -1;
    else
        vec[2] = atts.GetZMax();

    sel->SetStops(vec);
    vec[0] = atts.GetXIncr();
    vec[1] = atts.GetYIncr();
    vec[2] = atts.GetZIncr();
    sel->SetStrides(vec);
    selID = rv->GetDataRequest()->AddDataSelection(sel);

    if (rv->GetDataRequest()->MayRequireNodes())
        rv->GetDataRequest()->TurnNodeNumbersOn();

    return rv;
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::PreExecute
//
//  Purpose:
//      Called before Execute, which in turn calls ExecuteData.
//
//  Programmer: Hank Childs
//  Creation:   June 29, 2002
//
// ****************************************************************************

void
avtIndexSelectFilter::PreExecute(void)
{
    avtPluginDataTreeIterator::PreExecute();
    atLeastOneThreadSuccessfullyExecuted = false;
    if (!GetInput()->GetInfo().GetValidity().GetZonesPreserved())
    {
        if (atts.GetXIncr()!=1 || atts.GetYIncr()!=1 || atts.GetZIncr()!=1)
        {
            avtCallback::IssueWarning("The data was already modified "
                           "before the index select operator was applied."
                           "  It is only possible to do increments of 1.");
        }
    }
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::PostExecute
//
//  Purpose:
//      Called after Execute (which in turn called ExecuteData).
//
//  Programmer: Hank Childs
//  Creation:   June 29, 2002
//
// ****************************************************************************

void
avtIndexSelectFilter::PostExecute(void)
{
    avtPluginDataTreeIterator::PostExecute();
    
    int topoDim =
      GetInput()->GetInfo().GetAttributes().GetTopologicalDimension();

    int dstype = GetInput()->GetInfo().GetAttributes().GetMeshType();

    // If the initial mesh is unstructure then cheating and processing
    // an unstrucrutured mesh so skip checking for a topology
    // dimension reduction.
    if (topoDim > 0 && dstype != AVT_UNSTRUCTURED_MESH)
    {
      if (atLeastOneThreadSuccessfullyExecuted)
      {
        int newdim =
          GetInput()->GetInfo().GetAttributes().GetTopologicalDimension();

        switch (atts.GetDim())
        {
        case IndexSelectAttributes::ThreeD:
          if (atts.GetZMin() == atts.GetZMax())
            --newdim;
          
          // FALLTHRU
          
        case IndexSelectAttributes::TwoD:
          if (atts.GetYMin() == atts.GetYMax())
            --newdim;
          
          // FALLTHRU
          
        case IndexSelectAttributes::OneD:
          if (atts.GetXMin() == atts.GetXMax())
                --newdim;
          break;
              
        default:
          EXCEPTION0(ImproperUseException);
        }
        
        newdim = (newdim < 0 ? 0 : newdim);
        GetOutput()->GetInfo().GetAttributes().SetTopologicalDimension(newdim);
      }
    }

    if ( GetOutput()->GetInfo().GetAttributes().GetTopologicalDimension() > 0 &&
        (atts.GetXWrap() || atts.GetYWrap() || atts.GetZWrap()))
    {
        Replicate();

        avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
        outAtts.GetOriginalSpatialExtents()->Clear();
        outAtts.GetDesiredSpatialExtents()->Clear();
        outAtts.GetActualSpatialExtents()->Clear();

        double bounds[6];
        avtDataset_p ds = GetTypedOutput();
        avtDatasetExaminer::GetSpatialExtents(ds, bounds);
        outAtts.GetThisProcsOriginalSpatialExtents()->Set(bounds);
    }
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::ReleaseData
//
//  Purpose:
//      Releases the problem size data associated with this filter.
//
//  Programmer: Hank Childs
//  Creation:   September 10, 2002
//
//  Modifications:
//
//    Hank Childs, Fri Mar  4 08:12:25 PST 2005
//    Do not set outputs of filters to NULL, since this will prevent them
//    from re-executing correctly in DLB-mode.
//
//    Hank Childs, Fri Mar 11 07:37:05 PST 2005
//    Fix non-problem size leak introduced with last fix.
//
//    Kathleen Bonnell, Mon Jan 30 15:10:26 PST 2006
//    Handle vtkMaskPoints. 
//
// ****************************************************************************

void
avtIndexSelectFilter::ReleaseData(void)
{
    avtPluginDataTreeIterator::ReleaseData();
    lspace.clear();
}

// ****************************************************************************
//  Method: avtIndexSelectFilter::UpdateDataObjectInfo
//
//  Purpose:
//    Indicates that original nodes are required for Pick, and that
//    original zones cannot be used with Pick.
//
//  Programmer: Kathleen Bonnell 
//  Creation:   August 4, 2005 
//
//  Modifications:
//    Kathleen Bonnell, Mon May  1 08:57:41 PDT 2006
//    Changed call from OrigNodes to OrigElements, indicating that either
//    nodes or zones are required, or both. 
//
//    Hank Childs, Sat Mar  3 16:28:16 PST 2007
//    Put in about data attributes that we removed ghost data.
//
//    Kathleen Biagas, Tue Apr  1 14:17:22 PDT 2014
//    Invalidate zones.
//
// ****************************************************************************

void
avtIndexSelectFilter::UpdateDataObjectInfo(void)
{
    //
    // Node Pick returns wrong results on an Index selected plot unless it has 
    // original node numbers.  (The points are not transformed, but their 
    // numbering is probably different.  So set a flag that specifies that
    // they are needed for pick. 
    //
    GetOutput()->GetInfo().GetAttributes().SetOrigElementsRequiredForPick(true);

    //
    // Zone Pick CANNOT use Original zone numbers on an Index selected plot
    // because it may be the case that MANY original zones map to a SINGLE
    // current zone.  So set a flag specifying that the original zones
    // array CANNOT be used with pick.
    //
    GetOutput()->GetInfo().GetAttributes().SetCanUseOrigZones(false);
    GetOutput()->GetInfo().GetAttributes().SetContainsGhostZones(AVT_NO_GHOSTS);

    //
    // Indicate zone numbering probably changed
    //
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
}


// ****************************************************************************
//  Method:  avtIndexSelectFilter::FilterUnderstandsTransformedRectMesh
//
//  Purpose:
//    If this filter returns true, this means that it correctly deals
//    with rectilinear grids having an implied transform set in the
//    data attributes.  It can do this conditionally if desired.
//
//  Arguments:
//    none
//
//  Programmer:  Jeremy Meredith
//  Creation:    February 15, 2007
//
// ****************************************************************************

bool
avtIndexSelectFilter::FilterUnderstandsTransformedRectMesh()
{
    // Index select is based on logical coordinates only.
    return true;
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::VerifyInput
//
//  Purpose:
//    Throw an exception if user-selected domain is out of range. 
//
//  Programmer: Kathleen Bonnell 
//  Creation:   June 7, 2007 
//
//  Modifications:
//    Kathleen Bonnell, Thu Jun 21 16:31:59 PDT 2007
//    Determine amrLevel during set validation.
//
//    Eric Brugger, Wed Dec  3 08:23:58 PST 2008
//    I modified the routine to set groupCatergory to true when "Use Whole
//    Collection" was set, so that it would would index select based on the
//    whole mesh (using base_index) and not on a per block basis.
//
//    Hank Childs, Mon Dec 14 16:55:10 PST 2009
//    Updated for new SIL interface.
//
// ****************************************************************************

void
avtIndexSelectFilter::VerifyInput()
{
    if (atts.GetUseWholeCollection() || atts.GetSubsetName() == "Whole")
    {
        groupCategory = true;
        return;
    }

    std::string category = atts.GetCategoryName();
    std::string subset = atts.GetSubsetName();
    avtSILRestriction_p silr = GetOriginatingSource()->
        GetFullDataRequest()->GetRestriction();

    int setID, collectionID;
    TRY
    {
        collectionID = silr->GetCollectionIndex(category, silr->GetTopSet());
        setID = silr->GetSetIndex(subset, collectionID);
        avtSILCollection_p coll = silr->GetSILCollection(collectionID);

        if (coll->GetRole() != SIL_DOMAIN && coll->GetRole() != SIL_BLOCK)
        {
            //
            //  May occur if user types in a category name.
            //
            EXCEPTION1(InvalidCategoryException, category.c_str()); 
        }

        bool validSet = false;
        int nEls = coll->GetNumberOfSubsets();
        for (int i = 0; i < nEls && !validSet; i++)
        {
            validSet = (setID == coll->GetSubset(i));
            if (validSet && coll->GetRole() == SIL_BLOCK)
                amrLevel = i;
        }

        if (!validSet)
        {
            //
            //  May occur if user types in a set name.
            //
            EXCEPTION2(InvalidSetException, category.c_str(), subset.c_str());
        }

        if (coll->GetRole() == SIL_BLOCK)
        {
            groupCategory = true;
        }
    }
    CATCH(InvalidVariableException)
    {
        //
        //  SIL could not match category name or subset name to an id.
        //
        RETHROW; 
    }
    ENDTRY
}

vtkDataArray *avtIndexSelectFilter::GetCoordinates( vtkRectilinearGrid *grid,
                                                    unsigned int coor)
{
  if( coor == 0 )
    return grid->GetXCoordinates();
  else if( coor == 1 )
    return grid->GetYCoordinates();
  else if( coor == 2 )
    return grid->GetZCoordinates();
  else
    return 0;
}


void avtIndexSelectFilter::SetCoordinates( vtkRectilinearGrid *grid,
                                           vtkDataArray *coordinates,
                                           unsigned int coor)
{
  if( coor == 0 )
    grid->SetXCoordinates(coordinates);
  else if( coor == 1 )
    grid->SetYCoordinates(coordinates);
  else if( coor == 2 )
    grid->SetZCoordinates(coordinates);
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::Replicate
//
//  Purpose:
//    Replicates the first slice to the last slice for wrapping.
//
//  Programmer: Hank Childs/Allen Sanderson
//  Creation:   April 14, 2010
//
//  Modifications:
//    Kathleen Biagas, Tue Aug 21 16:14:59 MST 2012
//    Preserve coordinate type.
//
//    Kathleen biagas, Tue June 9 09:36:27 MST 2015
//    Changed signatu
// ****************************************************************************

void
avtIndexSelectFilter::Replicate()
{
    avtDataset_p ds = GetTypedOutput();

#ifdef PARALLEL
    int gd[6];
    MPI_Allreduce(&globalDims[0], &gd[0], 3, MPI_INT, MPI_MIN, VISIT_MPI_COMM);
    MPI_Allreduce(&globalDims[3], &gd[3], 3, MPI_INT, MPI_MAX, VISIT_MPI_COMM);
    for (int i = 0; i < 6; ++i)
        globalDims[i] = gd[i];
#endif
    int myRank = PAR_Rank();

    lspace.clear();
    int nblocks = 0;
    vtkDataSet **blocks = ds->GetDataTree()->GetAllLeaves(nblocks);

    vector<int> domainIds;
    vector<string> labels;
    ds->GetDataTree()->GetAllDomainIds(domainIds);
    ds->GetDataTree()->GetAllLabels(labels);
    bool wrap[3] = {atts.GetXWrap(), atts.GetYWrap(), atts.GetZWrap()};

    bool haveNewOutput = false;

    int dataObjectType = -1;

    bool okay = true;
    for (int n = 0; n < nblocks && okay; ++n)
    {
        vtkDataSet *ds = blocks[n]; 
        dataObjectType = ds->GetDataObjectType();
        int dims[3];
        vtkIntArray *bi = vtkIntArray::SafeDownCast(
                              ds->GetFieldData()->GetArray("base_index"));
        if (!bi)
        {
            okay = false;
            break;
        }
        int *mins = (int*)bi->GetVoidPointer(0);
        int max[3];
        if (dataObjectType == VTK_RECTILINEAR_GRID)
        {
            vtkRectilinearGrid *rgrid = vtkRectilinearGrid::SafeDownCast(ds);
            rgrid->GetDimensions(dims);
            bool doReplicate = false;
            for (int i = 0; i < 3; ++i)
            {
                max[i] = mins[i]+dims[i] -1;
                doReplicate |= (wrap[i] && max[i] == globalDims[i+3]);
            }
            if (doReplicate)
            {
                blocks[n] = Replicate(rgrid, wrap, dims, max);
                haveNewOutput = true;
            }
        }
        else if (dataObjectType == VTK_STRUCTURED_GRID)
        {
            vtkStructuredGrid::SafeDownCast(ds)->GetDimensions(dims);
            avtIndexSelectFilter::LogicalSpaces ls;
            for (int i = 0; i < 3; ++i)
            {
                max[i] = mins[i]+dims[i] -1;
            }
            ls.SetMins(mins);
            ls.SetMaxs(max);
            ls.rank = myRank;
            ls.block =  domainIds[n];
            lspace.push_back(ls);
        }
        else
        {
            okay = false;
            break;
        }

    }
#ifdef PARALLEL
    okay = (bool) UnifyMaximumValue((int)okay);
    dataObjectType = UnifyMaximumValue(dataObjectType);
#endif
    if (!okay)
    {
        lspace.clear();
        return;
    }

    if (dataObjectType == VTK_STRUCTURED_GRID)
    {
#ifdef PARALLEL
        int tags[2];
        GetUniqueMessageTags(tags, 2);
        int mpiNBTag   = tags[0];
        int mpiDataTag = tags[1];
        vector<avtIndexSelectFilter::LogicalSpaces> maxSpaces;
        vector<avtIndexSelectFilter::LogicalSpaces> minSpaces;
        if (myRank == 0)
        {
            MPI_Status stat;
            MPI_Status stat2;
            for (int i = 1; i < PAR_Size(); ++i)
            {
                int nb;
                MPI_Recv(&nb, 1, MPI_INT, MPI_ANY_SOURCE, mpiNBTag,
                         VISIT_MPI_COMM, &stat);
                for (int j = 0; j < nb; ++j)
                {
                    int data[8];
                    MPI_Recv(data, 8, MPI_INT, stat.MPI_SOURCE,
                             mpiDataTag, VISIT_MPI_COMM, &stat2);
                    avtIndexSelectFilter::LogicalSpaces ls;
                    ls.rank  = data[0];
                    ls.block = data[1];
                    ls.SetMins(&data[2]);
                    ls.SetMaxs(&data[5]);
                    lspace.push_back(ls);
                }
            }
        } 
        else
        {
            MPI_Send(&nblocks, 1, MPI_INT, 0, mpiNBTag, VISIT_MPI_COMM);
            for (int i = 0; i < nblocks; ++i)
            {
                MPI_Send(&(lspace[i].rank), 8, MPI_INT, 0, mpiDataTag, 
                           VISIT_MPI_COMM);
            } 
        } 
        
        Barrier();
#endif

        for (int i = 0; i < 3 && okay; ++i)
        {
            if(!wrap[i])
                continue;

        if (myRank == 0)
        {
            for (size_t j = 0; j < lspace.size(); ++j)
            {
                if (lspace[j].MatchesMaxAt(i, globalDims[i+3]))
                {
                    int findMin[3];
                    for (int k = 0; k < 3; ++k)
                    {
                        if (k == i)
                            findMin[k] = globalDims[i];
                        else
                            findMin[k] = lspace[j].mins[k];
                    }
                    for (size_t k = 0; k < lspace.size(); ++k)
                    {
                        if (lspace[k].MatchesMins(findMin))
                        {
#ifdef PARALLEL
                            maxSpaces.push_back(lspace[j]);
                            minSpaces.push_back(lspace[k]);
#else
                            int maxBlock = lspace[j].block;
                            int minBlock = lspace[k].block;
                            blocks[maxBlock] = Replicate(i, blocks[minBlock],
                                                         blocks[maxBlock]);
                            haveNewOutput = true;
#endif
                            break;
                        }
                    }
                }
            }
        } // PAR_Rank
#ifdef PARALLEL
        Barrier();
        int numMatches = (int)maxSpaces.size();
        BroadcastInt(numMatches);
        maxSpaces.resize(numMatches);
        minSpaces.resize(numMatches);
        for (int j = 0; j < numMatches; ++j)
        {
            BroadcastIntArray(&(maxSpaces[j].rank), 8);
            BroadcastIntArray(&(minSpaces[j].rank), 8);
        }
        for (size_t j = 0; j < maxSpaces.size(); ++j)
        {
            int tags[2];
            GetUniqueMessageTags(tags, 2);
            int mpiSizeTag = tags[0];
            int mpiDSTag   = tags[1];
            size_t maxBlock = 0;
            size_t minBlock = 0;
            if (myRank == maxSpaces[j].rank && myRank == minSpaces[j].rank)
            {
                // Both min and max blocks reside on this processor
                for (size_t k = 0; k < domainIds.size(); ++k)
                {
                    if (domainIds[k] == maxSpaces[j].block)
                        maxBlock = k;
                    if (domainIds[k] == minSpaces[j].block)
                        minBlock = k;
                }
                blocks[maxBlock] = Replicate(i, 
                    blocks[minBlock], 
                    blocks[maxBlock]);
                haveNewOutput = true;
            }
            else if (myRank == maxSpaces[j].rank)
            {
                // Need to recv the vtkDataSet from proc with min
                for (size_t k = 0; k < domainIds.size(); ++k)
                {
                    if (domainIds[k] == maxSpaces[j].block)
                        maxBlock = k;
                }
                MPI_Status stat;
                int size;
                vtkDataSetReader *reader = vtkDataSetReader::New();
                reader->ReadFromInputStringOn();
                MPI_Recv(&size, 1, MPI_INT, minSpaces[j].rank, mpiSizeTag,
                         VISIT_MPI_COMM, &stat);
                char *str = new char[size];
                MPI_Recv(str, size, MPI_CHAR, minSpaces[j].rank, mpiDSTag,
                          VISIT_MPI_COMM, &stat);
                vtkCharArray *charArray = vtkCharArray::New();
                charArray->SetArray((char*)str, size, 1);
                reader->SetInputArray(charArray);
                reader->Update();

                // Do the replication
                blocks[maxBlock] = Replicate(i, 
                    reader->GetOutput(), 
                    blocks[maxBlock]);
                haveNewOutput = true;

                // cleanup;
                delete [] str;
                reader->Delete();
                charArray->Delete();
              
            }
            else if (myRank == minSpaces[j].rank)
            {
                for (size_t k = 0; k < domainIds.size(); ++k)
                {
                    if (domainIds[k] == minSpaces[j].block)
                        minBlock = k;
                }
                // Need to send the vtkDataSet to proc with max
                vtkDataSetWriter *writer = vtkDataSetWriter::New();
                writer->WriteToOutputStringOn();
                writer->SetFileTypeToBinary();
                writer->SetInputData(blocks[minBlock]);
                writer->Write();
                int size = writer->GetOutputStringLength();
                char *str = writer->RegisterAndGetOutputString();
                MPI_Send(&size, 1, MPI_INT, maxSpaces[j].rank, mpiSizeTag,
                         VISIT_MPI_COMM);
                MPI_Send(str, size, MPI_CHAR, maxSpaces[j].rank, mpiDSTag,
                         VISIT_MPI_COMM);
                delete [] str;
                writer->Delete();
            }
        }
#endif
        } // wrap direction
    }


    if (haveNewOutput)
    {
        avtDataTree_p tree;
        if (!labels.empty())
        {
            tree = new avtDataTree(nblocks, blocks, domainIds, labels);
        }
        else 
        {
            tree = new avtDataTree(nblocks, blocks, domainIds);
        }
        SetOutputDataTree(tree);
    }

}


// ****************************************************************************
//  Method: avtIndexSelectFilter::Replicate
//
//  Purpose:
//    Replicates the first slice to the last slice for wrapping.
//
//  Arguments:
//    wrap      The wrap direction (0->i, 1->j, 2->k).
//    min_ds    The dataset with the minimum slice for the wrap direction.
//    max_ds    The dataset with the maximum slice for the wrap direction.
//
//  Programmer: Hank Childs/Allen Sanderson
//  Creation:   April 14, 2010
//
//  Modifications:
//    Kathleen Biagas, Tue Aug 21 16:14:59 MST 2012
//    Preserve coordinate type.
//
//    Kathleen biagas, Tue June 9 09:36:27 MST 2015
//    Changed signature (to handle multi-block and parallel), this method
//    now only handles vtkStructuredGrid.
//
// ****************************************************************************

vtkDataSet *
avtIndexSelectFilter::Replicate(int wrap, vtkDataSet *min_ds,
                                          vtkDataSet *max_ds)
{
    vtkStructuredGrid *min_sgrid = vtkStructuredGrid::SafeDownCast(min_ds);
    vtkStructuredGrid *max_sgrid = vtkStructuredGrid::SafeDownCast(max_ds);
    if (min_sgrid == NULL || max_sgrid == NULL)
    {
        debug3 << "IndexSelect::Replicate did not receive structured "
               << "grid input." << endl;
        return max_ds;
    }

    int   dims_in[3] = {0,0,0};
    int   dims_out[3] = {0,0,0};

    vtkDataSet *out_ds = max_ds;

    // Learn about the input grid
    max_sgrid->GetDimensions(dims_in);

    dims_out[0] = dims_in[0];
    dims_out[1] = dims_in[1];
    dims_out[2] = dims_in[2];
    dims_out[wrap]++;

    vtkPoints *pts = max_sgrid->GetPoints();
    vtkPoints *min_pts = vtkStructuredGrid::SafeDownCast(min_ds)->GetPoints();
    vtkPoints *new_pts = vtkPoints::New(pts->GetDataType());
    vtkStructuredGrid *out_sg = vtkStructuredGrid::New();
    out_sg->SetDimensions(dims_out);
    out_sg->SetPoints(new_pts);
    new_pts->Delete();
    new_pts->SetNumberOfPoints(dims_out[0]*dims_out[1]*dims_out[2]);

    // Copy the original points over.
    for (int i=0; i<dims_out[0]; ++i)
    {
        int i0 = i % dims_in[0];
    
        for (int j=0; j<dims_out[1]; ++j)
        {
            int j0 = j % dims_in[1];
     
            for (int k=0; k<dims_out[2]; ++k)
            {
                int k0 = k % dims_in[2];
                int idx_in  = k0*dims_in[1]*dims_in[0] + j0*dims_in[0] + i0;
                int idx_out = k*dims_out[1]*dims_out[0] + j*dims_out[0] + i;
                if (i < dims_in[0] && j < dims_in[1] && k < dims_in[2]) 
                    new_pts->SetPoint(idx_out, pts->GetPoint(idx_in));
                else
                    new_pts->SetPoint(idx_out, min_pts->GetPoint(idx_in));
            }
        }
    }

    if( out_ds != max_ds )
        out_ds->Delete();

    out_ds = out_sg;

    CopyData(max_ds, out_ds, dims_in, dims_out);

    return out_ds;
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::Replicate
//
//  Purpose:
//    Replicates the first slice to the last slice for wrapping.
//
//  Arguments:
//    rgrid     The input dataset.
//    wrap      Whether or not to wrap in i,j,k.
//    in_dims   The dimensions of the input grid.
//
//  Programmer: Hank Childs/Allen Sanderson
//  Creation:   April 14, 2010
//
//  Modifications:
//    Kathleen Biagas, Tue Aug 21 16:14:59 MST 2012
//    Preserve coordinate type.
//
//    Kathleen Biagas, Tue June 9 09:36:27 MST 2015
//    Extracted from original method to only handle rectilinear grid.
//
//    Kathleen Biagas, Wed June 10 10:45:43 MST 2015
//    Fixed logic so that wrapping for a given direction doesn't happend
//    if the input's max dims do not match the global max dims.
//
// ****************************************************************************

vtkDataSet *
avtIndexSelectFilter::Replicate(vtkRectilinearGrid *rgrid, bool wrap[3],
    int dims_in[3], int max[3])
{
    vtkRectilinearGrid *out = rgrid;
    bool haveCopied = false;
    int dims_out[3] = {dims_in[0], dims_in[1], dims_in[2] };

    // Do each dimension individually.
    for( unsigned int d=0; d<3; ++d )
    {
        if( !wrap[d] )
            continue;

        // does this dataset have the max dims?
        if (max[d] != globalDims[d+3])
            continue;

        if (!haveCopied)
        {
            out = rgrid->NewInstance();
            out->DeepCopy(rgrid);
            haveCopied = true;
        }
        dims_out[d]++;
        vtkDataArray *coor = GetCoordinates(out, d);

        double lastVal = coor->GetTuple1(dims_in[d]-1);
        double prevVal = coor->GetTuple1(dims_in[d]-2);
        coor->InsertNextTuple1(lastVal + (lastVal-prevVal));

    }

    if (haveCopied)
    {
        out->SetDimensions(dims_out);
        CopyData(rgrid, out, dims_in, dims_out);
    }

    return out;
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::CopyData
//
//  Purpose:
//    Copies Point and Cell Data from input to output datasets, utilizing
//    input and output structured dimensions.
//
//  Arguments:
//    in_ds     The input dataset.
//    out_ds    The output dataset.
//    in_dims   Dimensions for input dataset.
//    out_dims  Dimensions for output dataset.
//
//  Programmer: Kathleen Biagas 
//  Creation:   June 9, 2015
//
//  Modifications:
//
// ****************************************************************************

void
avtIndexSelectFilter::CopyData(vtkDataSet *in_ds, vtkDataSet *out_ds, int dims_in[3], int dims_out[3])
{
    // Copy over the point data.
    vtkPointData *inPD  = in_ds->GetPointData();
    vtkPointData *outPD = out_ds->GetPointData();
    outPD->CopyAllocate(inPD, dims_out[0]*dims_out[1]*dims_out[2]);

    for (int i1=0; i1<dims_out[0]; ++i1)
    {
        int i0 = i1 % dims_in[0];

        for (int j1=0; j1<dims_out[1]; ++j1)
        {
            int j0 = j1 % dims_in[1];

            for (int k1=0; k1<dims_out[2]; ++k1)
            {
                int k0 = k1 % dims_in[2];

                int idx_in  = k0*dims_in[1]*dims_in[0] + j0*dims_in[0] + i0;
                int idx_out = k1*dims_out[1]*dims_out[0] + j1*dims_out[0] + i1;

                outPD->CopyData(inPD, idx_in, idx_out);
            }
        }
    }

    // Copy over the cell data.
    int xdims = (dims_in[0]-1 < 1 ? 1 : dims_in[0]-1);
    int ydims = (dims_in[1]-1 < 1 ? 1 : dims_in[1]-1);
    int zdims = (dims_in[2]-1 < 1 ? 1 : dims_in[2]-1);

    int xdims_out = (dims_out[0]-1 < 1 ? 1 : dims_out[0]-1);
    int ydims_out = (dims_out[1]-1 < 1 ? 1 : dims_out[1]-1);
    int zdims_out = (dims_out[2]-1 < 1 ? 1 : dims_out[2]-1);

    vtkCellData *outCD = out_ds->GetCellData();
    vtkCellData *inCD  = in_ds->GetCellData();
    outCD->CopyAllocate(inCD, xdims_out*ydims_out*zdims_out);

    for (int i1=0; i1<xdims_out; ++i1)
    {
        int i0 = i1 % xdims;

        for (int j1=0; j1<ydims_out; ++j1)
        {
            int j0 = j1 % ydims;

            for (int k1=0; k1<zdims_out; ++k1)
            {
                int k0 = k1 % zdims;

                int idx_in = k0*(ydims*xdims) + j0*xdims + i0;
                int idx_out = k1*(ydims_out*xdims_out) + j1*xdims_out + i1;
                outCD->CopyData(inCD, idx_in, idx_out);
            }
        }
    }
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces constructor
//
//  Purpse:  Helper class for 'wrap' option.
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

avtIndexSelectFilter::LogicalSpaces::LogicalSpaces()
{
    mins[0] = mins[1] = mins[2] = -1;
    maxs[0] = maxs[1] = maxs[2] = -1;
}

// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces destructor
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

avtIndexSelectFilter::LogicalSpaces::~LogicalSpaces()
{
}

// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces::SetMins
//
//  Purpose:
//    Sets the minimum values.
//
//  Arguments:
//    _min     The minimum values.
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

void
avtIndexSelectFilter::LogicalSpaces::SetMins(int *_min)
{
    mins[0] = _min[0];
    mins[1] = _min[1];
    mins[2] = _min[2];
}


// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces::SetMaxs
//
//  Purpose:
//    Sets the maximum values.
//
//  Arguments:
//    _max     The maximum values.
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

void
avtIndexSelectFilter::LogicalSpaces::SetMaxs(int *_max)
{
    maxs[0] = _max[0];
    maxs[1] = _max[1];
    maxs[2] = _max[2];
}

// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces::HasMinAt
//
//  Returns:
//    True if the minimum value at index 'idx' matches the passed value.
//
//  Arguments:
//    idx     The index.
//    _min    The minimum value to be matched.
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

bool
avtIndexSelectFilter::LogicalSpaces::HasMinAt(int idx, int _min)
{
    return idx < 3 && idx >= 0 && mins[idx] == _min;
}

// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces::MatchesMin
//
//  Returns:
//    True if the mins values matches passed values.
//
//  Arguments:
//    i_min     The value to match at position 0.
//    j_min     The value to match at position 1.
//    k_min     The value to match at position 2.
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

bool
avtIndexSelectFilter::LogicalSpaces::MatchesMins(int _min[3])
{
    return mins[0] == _min[0] && mins[1] == _min[1] && mins[2] == _min[2];
}
 

// ****************************************************************************
//  Method: avtIndexSelectFilter::LogicalSpaces::HasMaxAt
//
//  Returns:
//    True if the maximum value at index 'idx' matches the passed value.
//
//  Arguments:
//    idx     The index.
//    _max    The maximum value to be matched.
//
//  Programmer: Kathleen Biagas
//  Creation:   June 9, 2015
//
// ****************************************************************************

bool
avtIndexSelectFilter::LogicalSpaces::MatchesMaxAt(int idx, int _max)
{
    return idx < 3 && idx >= 0 && maxs[idx] == _max;
}

